DEE Sampler Setup

DEV Sampler Setup

Sampler Examples

2D Example - well separated

Mixture of 3 multivariate normal densities all with diagonal variance structure and \(\sigma^2 = 10\). The true means are \((-20,20)\), \((20,-20)\), and \((0,0)\), and \(\sum_{j}n_j = 30\).

DEV sampler

without split-merge step
# n_iter = 5000
# use_cores = min(parallel::detectCores(), nreps)
# sim_results_DEV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = FALSE
#         )}
#       )
# 
# saveRDS(object = sim_results_DEV_no_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_noSM_6_Dec_2023.rds")
sim_results_DEV_no_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_noSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEV_no_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEV_no_sm[[rep]]$k, 
                        group_assign = sim_results_DEV_no_sm[[rep]]$group_assign)
  print(p1)
}
## 
##    1    2    3    4    5    6    7 
##    2   15 3965  856  131   24    7

## 
##    1    2    3    4    5 
##    4    7 4535  428   26

## 
##    1    3    4    5    6    7    8 
##    1 1694 2980  295   24    5    1

## 
##    1    3    4    5    6    7 
##    2 4447  478   61   10    2

## 
##    1    2    3    4    5    6    7 
##    5    1  459 3600  805  119   11

## 
##    1    2    3    4    5    6    7 
##    6    2 1632 3008  335   11    6

## 
##    1    2    3    4    5    6    7 
##    1    4 4094  818   70   11    2

## 
##    1    2    3    4    5    6    7 
##    4    4 2260 2310  373   44    5

## 
##    1    2    3    4    5    6 
##    4    5 3833 1079   76    3

## 
##    1    2    3    4    5    6 
##    5   44 3902  924  111   14

with split-merge step
# sim_results_DEV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = TRUE, sm_iter = 10
#         )}
#       )
# 
# saveRDS(object = sim_results_DEV_with_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_withSM_6_Dec_2023.rds")
sim_results_DEV_with_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_withSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEV_with_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEV_with_sm[[rep]]$k, 
                        group_assign = sim_results_DEV_with_sm[[rep]]$group_assign)
  print(p1)
}
## 
##    1    2    3    4    5    6    7 
##    2    5 3399 1243  309   41    1

## 
##    1    2    3    4    5    6 
##    5   11 4512  440   29    3

## 
##    1    2    3    4    5    6    7 
##    1    3 4302  603   78   10    3

## 
##    1    2    3    4    5    6    7 
##    2   55 3601 1117  196   26    3

## 
##    1    2    3    4    5    6    7 
##    5    1 3743 1108  127   15    1

## 
##    1    2    3    4    5    6    7    8 
##    6    2 3504 1207  248   27    5    1

## 
##    1    2    3    4    5    6 
##    1    4 3627 1205  152   11

## 
##    1    2    3    4    5    6 
##    4    4 3288 1521  158   25

## 
##    1    3    4    5    7 
##    4 4691  299    5    1

## 
##    1    2    3    4    5    6    7 
##    5   25 4074  765  120   10    1

DEE sampler

without split-merge step
# sim_results_DEE_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEE(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = FALSE
#         )}
#       )
# 
# saveRDS(object = sim_results_DEE_no_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_noSM_6_Dec_2023.rds")
sim_results_DEE_no_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_noSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEE_no_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEE_no_sm[[rep]]$k, 
                        group_assign = sim_results_DEE_no_sm[[rep]]$group_assign)
  print(p1)
}
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   1   7  23  94 230 450 669 710 712 610 496 351 280 172 111  49  17  14   3

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   2  14  71 216 449 620 740 715 682 518 333 262 176 115  50  20  12   2   1 
##  21 
##   1

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   1  11  67 261 538 724 739 666 593 392 301 234 170 138  80  46  26   8   4   1

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   2   1  13  91 196 359 574 686 756 650 516 457 311 191 101  66  21   7   2

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   6  43 144 327 523 627 734 640 611 494 331 261 131  67  41  14   4   1

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   3   3   1  20  92 277 515 723 805 680 626 419 325 206 147  90  49  12   6   1

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   1   5  42 146 387 620 845 800 722 488 338 245 158 102  56  24  14   3   1 
##  21 
##   2

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
##   4   3   2  25 126 376 621 804 787 698 522 387 280 159  95  49  36  22   4

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
##   4  61 217 445 677 780 781 728 538 333 207 118  77  22  12

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   5   3   9  32 109 259 468 680 765 715 574 476 347 231 135  85  53  39  10   5

with split-merge step
# sim_results_DEE_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEE(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = TRUE, sm_iter = 10
#         )}
#       )
# 
# saveRDS(object = sim_results_DEE_with_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_withSM_6_Dec_2023.rds")
sim_results_DEE_with_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_withSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEE_with_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k, 
                        group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
  print(p1)
}
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   1  14  30  87 273 468 725 750 663 580 492 365 221 161  96  47  17   6   3

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   2  78 111 259 432 632 690 748 652 441 351 281 135  81  66  23   8   7   2

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   1  24  69 250 536 781 788 672 539 406 272 231 161 118  67  46  22  12   3   2

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   2   4  41 164 371 600 755 706 660 494 410 275 226 141  84  38  19   6   1   2 
##  22 
##   1

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   4  39 151 336 541 682 689 675 553 458 358 222 144  83  32  18  12   2

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   3   3   4  26 114 300 541 716 820 688 560 425 309 209 136  76  45  17   7   1

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   1   4  29 180 403 718 855 755 627 498 320 236 167  84  67  35  13   3   4

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   4   3   7  48 189 413 647 805 762 673 500 355 244 149  99  55  25  13   8   1

## 
##   1   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
##   4  83 192 494 679 762 828 654 522 329 197 131  78  33   9   4   1

## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   5   3   4  43 137 345 556 700 741 658 514 430 303 201 148  96  62  30  19   5

UVV sampler

without split-merge step
sim_results_UVV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
      FUN = function(x){MVN_CRP_sampler_UVV(
        S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
        alpha = 1, r = 10, g = 1, h = 50, nu = 2,
        fix_r = FALSE, 
        mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
        lambda0 = diag(x = 15, nrow = 2),
        k_init = 1, diag_weights = FALSE, 
        verbose = FALSE, split_merge = FALSE
        )}
      )

saveRDS(object = sim_results_UVV_no_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_noSM_6_Dec_2023.rds")

# for(rep in 1:nreps){
#   print(table(sim_results_DEE_with_sm[[rep]]$k))
#   p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k, 
#                         group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
#   print(p1)
# }
with split-merge step
sim_results_UVV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
      FUN = function(x){MVN_CRP_sampler_UVV(
        S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
        alpha = 1, r = 10, g = 1, h = 50, nu = 2,
        fix_r = FALSE, 
        mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
        lambda0 = diag(x = 15, nrow = 2),
        k_init = 1, diag_weights = FALSE, 
        verbose = FALSE, split_merge = TRUE, sm_iter = 10
        )}
      )

saveRDS(object = sim_results_UVV_with_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_withSM_6_Dec_2023.rds")

2D Example - close together

Mixture of 3 multivariate normal densities all with diagonal variance structure and \(\sigma^2 = 10\). The true means are \((10,10)\), \((0,-5)\), , \((12,-2)\), and \((0,10)\), and \(\sum_{j}n_j = 50\).

DEV sampler

without split-merge step
# n_iter = 5000
# use_cores = min(parallel::detectCores(), nreps)
# sim_results_DEV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = FALSE
#         )}
#       )
# 
# saveRDS(object = sim_results_DEV_no_sm, file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

sim_results_DEV_no_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEV_no_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEV_no_sm[[rep]]$k, 
                        group_assign = sim_results_DEV_no_sm[[rep]]$group_assign)
  print(p1)
}
## 
##    1    2    3    4    5    6    7    8    9 
##   11  317 2558 1416  522  143   26    6    1

## 
##    1    2    3    4    5    6    7    8 
##   62  469 2370 1431  541  109   15    3

## 
##    1    2    3    4    5    6    7    8    9 
##  250  104 3296 1033  228   67   19    2    1

## 
##    1    2    3    4    5    6    7    8 
##  123  123 2814 1295  504  120   19    2

## 
##    1    2    3    4    5    6    7    8    9 
##  230  114 1643 1782  851  302   66   11    1

## 
##    1    2    3    4    5    6    7    8 
##  111  150 3301 1154  237   40    5    2

## 
##    1    2    3    4    5    6    7    8 
##   99   60 3836  795  152   50    7    1

## 
##    1    2    3    4    5    6    7    8    9 
##  359  604 1080 1566  897  375  100   18    1

## 
##    1    2    3    4    5    6 
##   51   20 4043  783  100    3

## 
##    1    2    3    4    5    6    7    8 
##  202  522 2331 1254  496  149   40    6

with split-merge step
# sim_results_DEV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = TRUE, sm_iter = 10
#         )}
#       )
# 
# saveRDS(object = sim_results_DEV_with_sm, 
#         file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")

sim_results_DEV_with_sm = readRDS(file = "../MCMC_Runs/conjDEVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEV_with_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEV_with_sm[[rep]]$k, 
                        group_assign = sim_results_DEV_with_sm[[rep]]$group_assign)
  print(p1)
}
## 
##    1    2    3    4    5    6    7    8 
## 1370  998 1134  920  401  145   26    6

## 
##    1    2    3    4    5    6    7 
## 2388 1425  795  275   95   19    3

## 
##    1    2    3    4    5    6    7 
## 2442 1132  789  457  147   30    3

## 
##    1    2    3    4    5    6    7    8 
## 1860 1193 1110  569  189   62   15    2

## 
##    1    2    3    4    5    6    7 
## 2857 1049  667  279  123   19    6

## 
##    1    2    3    4    5    6    7    8 
## 1515 1005 1300  730  320  105   19    6

## 
##    1    2    3    4    5    6    7 
## 3004 1250  502  174   49   18    3

## 
##    1    2    3    4    5    6    7    8 
## 2751 1083  629  355  125   40   14    3

## 
##    1    2    3    4    5    6    8 
## 1946 1168 1356  429   94    6    1

## 
##    1    2    3    4    5    6    7 
## 2193 1296  876  405  170   51    9

DEE sampler

without split-merge step
# sim_results_DEE_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEE(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = FALSE
#         )}
#       )
# 
# saveRDS(object = sim_results_DEE_no_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

sim_results_DEE_no_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEE_no_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEE_no_sm[[rep]]$k, 
                        group_assign = sim_results_DEE_no_sm[[rep]]$group_assign)
  print(p1)
}
## 
##    1    2    3    4    6    8    9   10   11   12   13   14   15   16   17   18 
##   15    4    1    1    1    2   21  115  343  729 1064 1016  787  486  258  103 
##   19   20   21   22   23 
##   36   13    3    1    1

## 
##   1   2   3   4   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##  29   2   1   2   1   3  18  72 225 540 756 952 944 706 431 209  68  32   7   2

## 
##   1   2   3   5   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22 
##  11   2   3   1   3  15  44 172 421 733 895 918 718 477 317 145  84  24  11   4 
##  23 
##   2

## 
##   1   2   3   4   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   2   2   1   1   1   1   8  49 196 508 831 980 917 650 419 227 114  59  24   9 
##  22 
##   1

## 
##    1    2    3    5    7    8    9   10   11   12   13   14   15   16   17   18 
##    8   10    1    1    1    2   14   89  358  706 1006 1036  773  518  271  130 
##   19   20   21   22 
##   53   16    5    2

## 
##    1    2    3    4    6    7    8    9   10   11   12   13   14   15   16   17 
##    2    1    1    1    1    1    7   43  213  465  832 1002  924  708  443  221 
##   18   19   20   21   22   23 
##   93   32    3    5    1    1

## 
##   1   2   3   4   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   2   2   1   1   2  11  42 174 409 692 830 914 758 511 318 181  89  40  16   4 
##  22  23 
##   2   1

## 
##   1   2   3   4   5   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 
##   1   3   2   1   1   1   1  22  72 253 556 852 921 887 633 439 210 105  35   5

## 
##    1    2    3    5    6    7    8    9   10   11   12   13   14   15   16   17 
##   20    1    2    2   15   53  189  467  814 1005  885  692  444  253  112   34 
##   18 
##   12

## 
##   1   2   3   4   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23 
##   8   1   1   1   1   6  28 133 363 660 909 992 807 541 333 145  49  19   1   1 
##  24 
##   1

with split-merge step
# sim_results_DEE_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_DEE(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50,
#         sigma_hyperprior = FALSE, fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         a = 1, b = 50,
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = TRUE, sm_iter = 10
#         )}
#       )
# 
# saveRDS(object = sim_results_DEE_with_sm, file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_withSM_6_Dec_2023.rds")

sim_results_DEE_with_sm = readRDS(file = "../MCMC_Runs/conjDEEsamp_minisimstudy_close_withSM_6_Dec_2023.rds")

for(rep in 1:nreps){
  print(table(sim_results_DEE_with_sm[[rep]]$k))
  p1 = make_k_traceplot(k = sim_results_DEE_with_sm[[rep]]$k, 
                        group_assign = sim_results_DEE_with_sm[[rep]]$group_assign)
  print(p1)
}

Some iterations failed.

UVV sampler

without split-merge step
# sim_results_UVV_no_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_UVV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50, nu = 2,
#         fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         lambda0 = diag(x = 15, nrow = 2),
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = FALSE
#         )}
#       )
# 
# saveRDS(object = sim_results_UVV_no_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

# sim_results_UVV_no_sm = readRDS(file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_noSM_6_Dec_2023.rds")

# for(rep in 1:nreps){
#   print(table(sim_results_UVV_no_sm[[rep]]$k))
#   p1 = make_k_traceplot(k = sim_results_UVV_no_sm[[rep]]$k, 
#                         group_assign = sim_results_UVV_no_sm[[rep]]$group_assign)
#   print(p1)
# }
with split-merge step
# sim_results_UVV_with_sm = parallel::mclapply(X = 1:10, mc.cores = use_cores,
#       FUN = function(x){MVN_CRP_sampler_UVV(
#         S = n_iter, seed = seeds[x], y = yreps[[x]]$y, 
#         alpha = 1, r = 10, g = 1, h = 50, nu = 2,
#         fix_r = FALSE, 
#         mu0 = matrix(round((colMeans(matrix(unlist(yreps[[x]]$y), ncol = 2))),0), ncol = 1), 
#         lambda0 = diag(x = 15, nrow = 2),
#         k_init = 1, diag_weights = FALSE, 
#         verbose = FALSE, split_merge = TRUE, sm_iter = 10
#         )}
#       )
# 
# saveRDS(object = sim_results_UVV_with_sm, file = "../MCMC_Runs/conjUVVsamp_minisimstudy_close_withSM_6_Dec_2023.rds")

This took too long…over 1.5 hours. Had to kill job for time. Need to reexamine and determine how to handle in future.